function [p,e,t,n_t,n_p,x_t,y_t,area]=mesh_refine(g)

[p,e,t]=initmesh(g);

ans=input('Need a Global Refinement? 1-yes,0-no=');
while ans==1
 [p,e,t]=refinemesh(g,p,e,t);

ans=input('Need Another Global Refinement? 1-yes,0-no=');
 if ans==0
  break;
 end
end %while

ans=input('Need Subdomain Refinement? 1-yes,0-no=');
while ans==1
 it=input('Input the subdomain number it=');
 [p,e,t]=refinemesh(g,p,e,t,it);

ans=input('Need another subdomain Refinement? 1-yes,0-no=');
 if ans==0
  break;
 end
end %while

n_t=size(t,2);  %n_t=total number of triangles, n_p - points
n_p=size(p,2);

%generate midpoint of each trangle
x_t=(p(1,t(1,:))+p(1,t(2,:))+p(1,t(3,:)))/3.;
y_t=(p(2,t(1,:))+p(2,t(2,:))+p(2,t(3,:)))/3.;

ans=input('Need a Selective Refinement? 1-yes,0-no=');
while ans==1
disp('Input the center (x_c y_c) and axes (r_x r_y) of the refining ellipse');
x_c=input('x_c=');  %center of an ellipse
y_c=input('y_c=');
r_x=input('r_x=');  %axes of the ellipse
r_y=input('r_y=');

j=1;

for i=1:n_p
   dist=(((x_t(i)-x_c)/r_x)^2+((y_t(i)-y_c)/r_y)^2);
if dist <= 1.00001
   id(j)=i;    %identify the elements inside the ellipse
   j=j+1;
end
end

[p,e,t]=refinemesh(g,p,e,t,id');

n_t=size(t,2);  %n_t=total number of triangles, n_p - points
n_p=size(p,2);

%----------generate midpoint of each trangle again---------
x_t=(p(1,t(1,:))+p(1,t(2,:))+p(1,t(3,:)))/3.;
y_t=(p(2,t(1,:))+p(2,t(2,:))+p(2,t(3,:)))/3.;

ans=input('Need Another Selective Refinement? 1-yes,0-no=');
 if ans==0
    break;
 end
end  %while
%*************************************************************

[area,a1,a2,a3]=pdetrg(p,t); %Generate area for each triangle  

ans=input('Want to plot the mesh? 1-yes,0-no=');
if ans==1
pdemesh(p,e,t);
end
